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Abstract 

We construct a program to calculate Feynman amplitudes at finite temper- 
ature in the real time Keldysh formalism using the symbolic manipulation 
program Mathematica. As an example, the usefulness of this program is 
demonstrated by proving the finite temperature Ward identity for QED in a 
second order effective theory. 
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I. INTRODUCTION 



There are two different methods for calculating Feynman amplitudes at finite temper- 
ature. Each of these methods presents technical difficulties The imaginary time 
formalism involves the calculation of green functions with imaginary time arguments. At 
the end of the calculation, one must perform an analytic continuation to obtain real time 
green functions. For higher n-point functions, these analytic continuations become increas- 
ingly difficult. In spite of this difficulty, the imaginary time formalism has been traditionally 
the most popular. 

The difficulty associated with the real time formalism is the doubling of degrees of free- 
dom. The real time integration contour involves two branches, one running from minus 
infinity to positive infinity just above the real axis, and one running back from positive 
infinity to negative infinity just below the real axis 0,^. All fields can take values on either 
branch of the contour and thus there is a doubling of the number of degrees of freedom. 
For example, the two point function becomes a two by two matrix. The four components of 
this matrix represent the four possible contractions of two field operators each of which can 
take values on either of the two branches of the real time contour. The physical two point 
functions (for example, the retarded and advanced two-point functions) can be extracted 
by taking apropriate combinations of the four components of this matrix. The procedure 
is similar for higher n-point functions. It is straightforward to show that this doubling of 
degrees of freedom is necessary to obtain finite green functions. 

The problem with the real time formalism is that the extra degrees of freedom become 
increasingly cumbersome to handle for higher n-point functions. Each line within a given 
Feynman diagram contains a Keldysh index at each end which takes values {1,2}, corre- 
sponding to the two branches of the closed time path contour. Indices that do not correspond 
to external legs are called internal indices, and must be summed over. The indices that cor- 
respond to external legs are called external indices. As a consequence of these indices, an 
n-point function has 2" components. These components obey one constraint equation, which 
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reduces the number of independent components to 2" — 1. In equilibrium, the KMS con- 
ditions impose additional constraints, reducing the number of independent components to 
2"-i - 1. 

In summary then, the imaginary time formalism has the advantage that the initial cal- 
culation is easier, but the disadvantage that analytic continuations must be done, which 
makes it difficult to extract the physical green functions. In contrast, in the real time for- 
malism, the initial calculation is more difficult, but there is a simple and natural procedure 
for extracting the physical green functions. In addition, the real time formalism can be gen- 
eralized to non-equilibrium situations. In this paper we describe a Mathematica calculation 
that does sums over internal indices, and takes physical combinations of external indices. 
Using this program substantially reduces the technical difficulties associated with the real 
time formalism, and makes it possible for us to exploit its advantages. 

This paper is organized as follows: In the section II we discuss the definitions of real time 
finite temperature green functions within the Keldysh formalism, and define our notation. 
In the section III, we describe the Mathematica program that we have written to do sum- 
mations over internal Keldysh indices, and take physical combinations of external indices. 
In the section III, we demonstrate the usefulness of this program by using it to prove the 
Ward identity for QED in a second order hard thermal loop (HTL) effective theory. In the 
last section, we give some conclusions. 

The program and the manuscript of this paper are available on the Los Alamos data 
base (|hep-ph / 9908438|) and at |http : / / theory. uwinnipeg. ca /users /EHI 



II. REAL TIME GREEN FUNCTIONS 

A. The two-point function 

We first consider the propagator. In real time, the propagator has 2^ = 4 components, 
since each of the two fields can take values on either branch of the contour. Thus, the 
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propagator can be written as a 2 x 2 matrix of the form 

], (1) 

V^21 D22J 

where Du is the propagator for fields moving along Ci, D12 is the propagator for fields 
moving from Ci to C2, etc. The four components are given by 

L'n(a:-y) = -^(T(0(x)0(|/))), 
Di2{x -y)^ -i{(t){y)(t){x)) , 
D2i{x-y) = -i{(j){x)(j){y)) , 

D22{x-y)^-i{f{cf>{x)cf>{y))), (2) 

where T is the usual time ordering operator, and T is the anti-chronological time ordering 
operator. These four components satisfy, 

Dn - D12 - D21 + L>22 = (3) 

as a consequence of the identity 9{x) + 0{—x) — 1. 

It is more useful to write the propagator in terms of the three functions 

Dr = Dn - D,2 , 
Da = Dn - D21 , 

Dp = Dn + D22 . (4) 
Dr and Da are the usual retarded and advanced propagators, satisfying 

Dn{x -y)- Da{x - y) = 0(y)]) , (5) 

and Dp is the symmetric combination 

DF{x-y) = -t{{<P{x),<P{y)}). (6) 
In momentum space these three propagators satisfy the KMS condition. 
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Df{P) = (1 + 2n{po))iDn{P) - Da{P)), (7) 

where 

(po ± ^ej — — 
and ri(po) is the Bose-Einstein distribution function, 

Equations (|^), @ are inverted by 

D22 = \{Df-Da-Dr). (10) 

These equations can be written in a more convenient notation as 

where the outer product of the column vectors is to be taken. 

Similar relations can be obtained for the IPI two-point function, or the polarization 
tensor, which is obtained by amputating the external legs from the propagator. The Dyson 
equation gives 

iD{p) = tDoip) + tDoip) (-zn(p)) tD{p) . (12) 
The analogues of (0) and (H) are 





= nn 


+ ni2 


Ha 


= nn 


+ n2i 




= nn 


+ n22 



(13) 



and 



Hii + + U21 + = . (14) 
The analogues of ([Tl|) and (0) are 

2 Uip) = Unip) Q + UM Q + (_\) (_\) , (15) 

= (1 + 2n(po)) (n^(p) - n^(p)) . (16) 



B. Three-Point Function 

In the real time formalism, the three-point function has 2'^ = 8 components. We denote 
the connected functions by F^^^ where {a, b,c = 1, 2}. In analogy to they are given by 
the following expressions [Q: 
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^122 
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^212 
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^221 
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y, 


z) 
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^222 
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= {f{<P{x)<Piy)<P{z))) . 



It is immediately obvious that one of these components is dependent on the others because 
of the identity 

EEE(-i)"^^^^-^rl = o (18) 

a=l b=l c=l 

which follows in the same way as (|^) from 6{x) + 6{—x) = 1. The seven combinations that 
we use are defined as 
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(19) 

In coordinate space we always label the first leg of the three-point function by x and call it 
the "incoming leg (i)", the third leg we label by z and call it the "outgoing leg (o)", and 
the second (middle) leg we label by y. Inserting the definitions (|T7| ) into (|19D one finds 
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3}, 02}) , 



(20) 

where we have used the obvious shorthands 0i = 0(x), 02 = (p{y), 03 = (p{z), and 612 = 
6{xq — yo), etc. The first three are the retarded three-point functions; F^^ is retarded with 
respect to Xq, F^^ is retarded with respect to Zq, and F§ is retarded with respect to y^. 

The IPI vertex functions are obtained from the connected functions by truncating ex- 
ternal legs. We will denote IPI vertex functions by F. We can write F as a tensor of the 
form 

-(:)(:)(:) 
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where the outer product of the column vectors is to be taken. For the IPI functions the 
analogues of (|T8|) and (|19]) are, 

2 2 2 

E E E Taftc = x + y + u + v + w + z = 0, (22) 

a=l 6=1 c=l 

and, 
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The IPI vertex functions r(Pi,P2,-P3) are related to the connected vertex functions 
r^(Pi,P2,P3) as follows: 
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where we have used the notation Dr{Pi) = ri, Df{P2) = /2, etc. 

For calculational purposes, we want to obtain a decomposition of the IPI three-point 
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function in terms of the seven functions p3|), in analogy to ([TT|) for the two-point function. 
Inverting (^) we obtain 

--(;)(■,)(;)".(-.)(;)(;)".(;)(;)(■,) 
--(■,)(:)('.)"-(;)(',)(■.)""(-',)(-',)(;)-(-',)(■,)(■.) 



C. Four-Point Function 

The connected four point function is given by the contour ordered expectation value, 

The IPI four-point function is obtained by truncating external legs and forms a 16 compo- 
nent tensor which we can write as the outer product of four two component vectors, 



M -- 

The retarded IPI four-point functions are given by 



Mm = Mini + Miu2 + Mu2i + M1211 + Mnaa + M1212 + M1221 + M1222 

= ^(x - y){u + v){w + z){s + t) 
Mr2 = Mini + Mni2 + Mnai + M2111 + Mn22 + ^2112 + M2121 + M2122 

= ^{x + y){u-v){w + z){s + t) (25) 

Mrs = Mini + Min2 + M2111 + M1211 + M2112 + M1212 + M2211 + M2212 
1 

= -(x + y){u + v){w - z){s + t) 

Mm = Mini + M2111 + Mn2i + M1211 + M2121 + M2211 + M1221 + M2221 
= ^(x + y){u + v){w + z){s - t) 

where we have used the relation 

2 

^ Mal>cd = Q. (26) 
a,f),c,d=l 
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The other combinations we will define as, 
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We use the decomposition of the four point vertex: 




D. Rules for calculating Feynman Amplitudes 

The rules for handeling Keldysh indices are as follows: (for details see Ref. p): 
1) Bare vertices carry a factor 
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2) Internal indices are to be summed over. In terms of the column vectors ocurring in (^Tj), 
(P^) and (p8D this means that one adds the product of the upper components to the product 
of the lower components of all column vectors carrying the same internal index. The product 
of any number of vectors carrying the same internal index gives a scalar: 

= X1X3 + X2X4 . (30) 

3) For external indices the product of any number of column vectors carrying the same 
index is defined to be another column vector whose upper (lower) component is given by 
the product of upper (lower) components of the original vectors: 

xi \ fx3\ _ /xixs 

X2J \X4j \X2X4 





III. THE PROGRAM 

The program was written entirely in Mathematica 3.0. Mathematica was chosen be- 
cause of its powerful numeric algorithms and ability to perform operations on sets The 
program performs contractions on Keldysh indices. It works for diagrams with three or four 
point interactions, with up to four external legs. Any number of loops can be considered, 
and any number of the vertices can be corrected vertices. Corrected vertices may be neces- 
sary when using an effective theory that involves a reorganized perturbation theory which is 
obtained by a resummation. For example, in the hard thermal loop approximation, it is not 
sufficient to consider only bare vertices. All fields are treated as scalars and the coefficient 
is calculated by assuming that a bare three-point vertex carries a factor —ig, and a bare 
four-point vertex carries a factor — zA. Sign conventions for the IPI functions are as shown 
in Fig. [1] . When non-scalar fields are involved, additional factors (such as traces over Dirac 
matrices for fermions, or contractions of projection operators for photons in a given gauge) 
must be calculated by hand. Also note that each line carries a factor of i which means 
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that if gauge boson propagators are present, usual conventions require the insertion of an 
additional minus sign if there is an odd number of gauge boson propagators. In addition, 
the sign for the gauge boson polarization tensor must be changed (see Fig. [1]). The user of 
the program must assign momenta to propagators and external legs so that the conservation 
of momentum is satisfied. Keldysh indices must also be assigned. When the program is 
executed, a number of input parameters are requested. We describe below the data entry 
process using the example shown in Fig. [2] . 

1) . The number of external legs is entered (3). For each external leg, the momentum and 
index are recorded in the following form, 

{P,a}; {K,b}; {PK,c} . 

Note that the momentum is entered without signs or spaces. The only purpose of these 
variables is to remind the user of the order he has chosen for the external legs. For the 
example used in this section, the program will calculate r(P, K, —P — K) since the order in 
which the external indices have been entered corresponds to this ordering of the momentum 
variables. 

2) . The number of internal indices (4) and a set containing these indices {d,e,f,g} is 
entered. 

3) . A set containing the loop variables is entered {R}. This set contains the independent 
momenta, and is used when the user wishes to have terms removed which are zero by contour 
integration. This point will be discussed further in section III-B. 

4) . The number of corrected three-point vertices is requested (1), and two further input 
boxes request the data pertaining to each corrected vertex. The first data group consists 
of a subscript, and the indices for the vertex {l,e, a, /}. The second group of data to be 
entered is the momenta corresponding to the previous indices {R, P, PR}. As before, these 
momenta are simply used as a record keeping device. To avoid clutter, the result is given 
without displaying explicitly the momentum dependence of the vertices. To remind the user 
of the choices he has made, the momentum dependence for each vertex is printed out. When 
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there are more than one corrected three-point vertices, the subscripts distinguish them. 

5) . The number of corrected four-point vertices (1) and their parameters are entered. The 
format is the same as for the corrected three-point vertices: {1, d, b, c, g}, {R, K, PK, PR}. 

6) . The number of indices that have a r associated with them is entered (0). If this number 
were non-zero, a further input box would request a set of these indices. 

7) . All of the data for each propagator is entered: the number of propagators (2), and 
each propagator's momentum, and initial and final index. In this case, the momenta are 
not merely used for book keeping purposes (as in the case of the external legs and corrected 
vertices) . The propagator momenta are used when eliminating terms that are zero by contour 
integration, and must be entered in a specific form. Each component of a propagator's 
momentum is entered seperately as a series of variables within a set bracket, and signs 
are included. If there is only one term in the momentum, it must still be enclosed by set 
brackets: 

{{R},d,e} 
{{P,R}J,g} 

8) . The user is now asked if he would hke the coefficient evaluated. To calculate the 
coefficient, the number of bare three-point vertices (0) and bare four-point vertices (0) must 
be entered. For three- and four-point functions, the coefficient is calculated using the formula 

coefficient = {^J {-igfi-iXy 

where p is the number of propagators, t is the number of bare three-point vertices, / is the 
number of bare four-point vertices, F is the number of corrected three-point vertices, M 
is the number of corrected four-point vertices, and in is the total number of indices. For 
two-point functions, an extra factor of —i is included (see Fig. [1]). 

9) . At this point, the user is asked which combination of external indices he would like 
evaluated. For example, we will choose to evaluate the retarded combination T^^P, K, —P — 
K). 
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10) . At the user's discretion, the result of the calculation can be displayed with or without 
the terms removed which are zero by contour integration. If the program is asked to remove 
these terms, it will remove terms in which, for any loop momentum, all poles of the propa- 
gators are on the same side of the real axis. Since the complex integral can be evaluated by 
choosing a semi-circle in the upper or lower half plane, these terms normally give zero. In 
some cases however (diagrams with tadpoles) , there are terms with all poles are on the same 
side of the real axis which do not give zero (because of contributions from the semi-circle 
at infinity). In addition, the program does not consider poles within the corrected vertices. 
When corrected vertices are present, or if there are tadpole pieces to the diagram, the user 
must tell the program not to remove any of the terms that it thinks will be zero. 

11) . The user can then choose to evaluate any other combination of external indices using 
the initial data. For example, he could choose to calculate Tpi^P, K, —P — K). 

The result is displayed in one of two ways. If there are less than four propagators, the 
result is simply shown as an unfactored polynomial. If there are four or more propagators, 
terms involving the first two propagators entered are factored out of the entire polynomial. 
The Mathematica operation Simplify[. . .] is used on the remaining terms. In addition, 
certain elements of the initial input are echoed as output to help the user detect typing 
mistakes. 

For the example shown in Fig. [1] the result is; 



where we have used the notation dR = d'^r/ (2n)'^, = Dr{R), fp^^ = Dp{P + R) etc. This 
example contains two propagators (each of which contains three terms) , one corrected three- 
point function (which contains seven terms), and one corrected four-point function (which 
contains 15 terms). To do the calculation by hand, one would have to evaluate 3^- 7 - 15 = 945 
terms. In fact, only five terms are non-zero. The program identifies these non-zero terms. 



Tr{P,R,-P-R) 




M := M(P + R,-P- K, K, -R) F := F(P, R, -P - R) 
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The use of a program of this type makes real time finite temperature calculations practical, 
and allows us to exploit the advantages of the real time formalism, one of which is the fact 
that one obtains physical green functions directly. 



IV. THE WARD IDENTITY FOR QED IN A SECOND ORDER HTL 

EFFECTIVE THEORY 

In a gauge theory, the Ward identities are a reflection of the gauge symmetry: if the 
theory is invarient under gauge transformations, then the green functions of the theory obey 
the Ward identities. It is well known that the QED Ward identities hold in a first order 
HTL effective theory. This result is a consequence of the fact that the HTL theory respects 
gauge invarience. In this section we will verify that the Ward identity for the three-point 
function and the electron self-energy is obeyed in a second order HTL effective theory. We 
will show that the Ward identity 

ir^r^ = -^e[S^(P) - S^(P + K)] (32) 

is satisfied by the diagrams shown in Fig. [3] and Fig. [4]. In these diagrams the solid dots 
are corrected vertices, which are obtained by adding the HTL vertex to the bare vertex, 
and the dotted lines are HTL propagators. We will work in real time and use the program 
described in section HI to perform the summations over Keldysh indices. This calculation 
would be prohibitively tedious using standard calculational techniques. 
We use the Coulomb gauge. The gauge boson propagator is given by, 

1 



00 



P - H 



00 

where we use the notation K^j, = {kQ,k), and Hqo and Ha are components of the HTL 
polarization tensor. For simplicity, we will consider only longtitudinal modes. To avoid the 
introduction of more notation, we will not explicitly distinguish these propagators from the 
bare propagators in section II. For example, we write 
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and 

/p = <[^p-«p]; < = l + 2n(po) (33) 

where n{po) is the Bose-Einstein distribution function @. Fermion propagators are written 
with tildes. For example, 

S^\P) = r-^R{P); r, = Sn{P) 

where S(P) is the hard thermal loop fermion self-energy. The symmetric propagator is 
defined as 0, 

fp = NF{P)P{rp-dp) (34) 
where Np{P) is constructed from the Fermi-Dirac distribution function, 

iV^^ = 1 - 2nf{po) ; n^ipo) = ^^^^ (35) 

To simplify notation, we suppress the Lorentz index on all vertices. For example, for the 
three-point vertex we write Fq := F. For the four-point vertex we write M^q := M^. In 
addition, all Dirac indices are suppressed. 

We will need to write four- and three-point functions with various momentum depen- 
dencies. Suppressing all indices for the moment, we make the following definitions, 

M(*) := M{P + R,-P- K, K, -R) 
M(**) := M(P, -P - K - R,K, R) 

._ Y{P,R,-P - R) 
YiB) ._ Y(p + K + R, -R, -K - P) (36) 
Yic) ._ p(p ^ ji x,-P - K - R) 

We will also need vertices with the signs of the momenta reversed. We use the following 
notation, 
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if r(^) = r(Pl,P2,P3) then T^^") = r(-Pl, -P2, -P3) (37) 

The four- and three-point HTL corrected vertices are related through the Ward identities. 
The identities that we will need are HI, 
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(38) 



Using the Mathematica program described in section III, we calculate the contribution 
from the diagrams in Fig. [3A] and Fig. [3B]: 

Tt\P, K,-P-K) = ydR[ (Mg)'^(r,Sp+,r^, + r J.+.F^ + /.S.+.F^) 

+ (M«)'^r,f,+pF^ + {M^2^rardr+pTt ] (39) 

Ff )(P, K,-P-K) = ^JdR[ (mJ^J )^(a,f,+fc+,F|„ + a J.+.+.F^ + /.f,+,+,F|J 

+ iMi}'^rar~ar^,^,T^ + (Mf ^r^f.+fc+.F^, ] (40) 

Using the Ward identity (|38|) we obtain, 

K^{Tt^+Tf^) = ydR{a + P) (41) 



where 



a — F^j (r^ap+^rFi + ''^rfp+r^R + frO'p+r^^o) + ^F f^r^r+p^R + ^ Fo '^r-C'r+p^Ro 

~[^Ro {(^T-f^'p+k+r^Fo + (^rfp+k+r^R + fr^p+k+r^Ri) (42) 

"t" i pi U.rar+p+fci Ji i- Pi I r' r+p+k^ Ril 
(3 = —Tj^-(rrdp^r^pj^ + Trfp+r^R + /rSp+rF^^) — T p^a-rdp-^-rT 

+^Ro{(^rfp+k+r^Fo + (^rfp+k+r^R + frrp+k+r^Ri) + ^Fi^k^p+k+r^Ri (43) 
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The Mathematica program gives the retarded and advanced self-energies in Fig. [4] as, 

S-R(-P ~^ ~ 2 I ^ \<^rrr+p+k^Fo + (^rfr+p+k^R + frrr+p+k^Ri) 

' Fi 1" r+p+fc-L Ri ' ^ F "r^r+p+k^ R J 

+ T^t\rfp+rTi + Tf~^ar~a,+rTl ] (44) 



Comparing (^21) and (0) we find that the contributions from the a term alone give, 

[K^{vf^ + rf only = -^e(S^(P) - S^(P + ir)) (45) 

Next we will show that the (3 terms ( ^31) cancel with the contributions from the diagram 
in Fig. [3C]. Using the Mathematica program as before, we obtain the contribution from 
Fig. [3C]. The result is as follows: 

Tf\P,K,-P-K)='-j dR 

[ I rU-p+r' p+k+T-yi- R jpL'-Fi'-Ri ~r ' rj p+r' p+k+ryl- R Jfi^-R >- Ri 
-t-arap+r^p+fc+rli ^ j^i RqI Fo + ^r^p+r/p+fc+r U _R Jm-*- Ro R l^^J 
^JrU-p+r' p+k+ry- R ) {m'- Ro Ri ^ ' r' p+r' p+k+r\'- Fo ) p.'- R Ri 
+Q'rCtp+rCtp+A:+r-U j^i /jo i J 

Contracting with K^^ and using the Ward identities for the HTL corrected three point 
vertices: 

^'^(^Fi^)lJ. ~ ^^-^p+k+ri^p+k+r ~ ^p+k+r) 
^^(^R = ^^^p+ri'^p+r " ^p+r) 

(47) 

we find that K^T'^^'* cancels exactly with (^Sp. As a result, we obtain from (P5|), 

ir^(r^") + rg) + rg)), = -.e(s^(P) - s^(p + ir)) (48) 



This is the Ward identity with full propagators and corrected vertices in QED at finite 
temperature. 

V. CONCLUSION 

We have constructed a program with Mathematica to evaluate Feynman amplitudes 
in the Keldysh formalism of real time finite temperature field theory. This formalism has 
recently gained increasing popularity because it avoids the need for analytical continuations 
that plagues the imaginary formalism, and it allows for a generalization to non-equilibrium 
situations. However, because of the extra degrees of freedom, calculations in the real time 
formalism can be extremely tedious, especially for higher n-point functions. We have written 
a Mathematica program that performs sums over Keldysh indices in the real time formal- 
ism. The program calculates physical Feynman amplitudes for any diagram with up to four 
external legs, with an arbitrary number of loops. Generalization to diagrams with more ex- 
ternal legs is straightforward. This program makes real time finite temperature calculations 
feasible for diagrams with complicated structure. In section III a diagram with 945 terms is 
calculated as an example. The program performs the summations over Keldysh indices and 
produces a result in which only five terms are non-zero. 

In order to demonstate the usefulness of this program, we have used it to prove the 
finite temperature QED Ward identity in a second order HTL effective theory. The relevant 
diagrams include corrected three- and four-point vertices and full propagators. To calculate 
them by hand would be extremely time consuming. Using the program developed in this 
paper, we are able to calculate these diagrams and verify the QED Ward identity with a 
minimum of effort. 
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Figure Captions 

Fig. 1. Conventions for the definitions of the vertex funtions. 

Fig. 2. Example of a three-point function used to describe input parameters. 

Fig. 3. Three-point functions that contribute to the Ward identity. 

Fig. 4. Two-point function that contributes to the Ward identity. 
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